/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */

package MF1DSolver;

import FEMShapeFunction1D.FEMShapeFunction1D1;
import FEMSolver1D.MainMatrixGenerator1;
import FEMSolver1D.MatrixGenerator;
import nonequidistant.FEMSolver1D.MatrixGenerator1;
import HashMatrix2.HashMatrix2;
import HashMatrix2.LUFactorization;
import Integration.GaussQuadrature;
import Integration.IFunction1D;
import MathLib.HashMatrix.Gauss;
import MathLib.HashMatrix.HashMatrix;
import MeshlessShapeFunctions.IKernelFunction;
import MeshlessShapeFunctions.JCubicSpline;
import MeshlessShapeFunctions.JGaussian;
import MeshlessShapeFunctions.JMeshlessInterpolationFunction1D;
import MeshlessShapeFunctions.JMeshlessInterpolationFunctionLiu;
import MeshlessShapeFunctions.JStKernel;
import MeshlessShapeFunctions.Particle;
import java.util.ArrayList;
import java.util.Iterator;
import java.util.logging.Level;
import java.util.logging.Logger;

import nonequidistant.FEMShapeFunctions1D.FEMShapeFunction1D1Pow;
import nonequidistant.FEMSolver1D.FEMResult;
import nonequidistant.FEMSolver1D.StatMatrixGenerator;

import org.jfree.data.xy.XYDataItem;
import org.jfree.data.xy.XYSeries;
import org.jfree.data.xy.XYSeriesCollection;

import syntaxanalyser.PostfiksX;

import jmeshless.NewOkCancelDialog;


/**
 *
 * @author mark_o
 */
public class MF1SolverLagrangeCoupled {
	public class optimizationResult{
		protected int N; // meshFree
		protected int M; // FEM
		protected double divPoint;
		
		public optimizationResult(int N, int M, double divPoint) {
			this.divPoint = divPoint;
			this.M = M;
			this.N = N;
		}
		
		public double getDivPoint() {
			return divPoint;
		}
		
		public int getM() {
			return M;
		}
		
		public int getN() {
			return N;
		}
		
		
	}


	public static double dirihletB = 0;
	public static double dirihletA = 0;
	public static String f="100";
	public static String p="1";
	public static String q="0";
	public static String w="100";
	
	
    //результуючі матриця та вектор
    HashMatrix A;
    double[] b;

    //фпроксимація задачі
    HashMatrix2 Afem;
    HashMatrix2 Amf;
    //матриці мас
    ArrayList<Double> Bfem;
    ArrayList<Double> Bmf;
    ArrayList<Double> BPoxfem;
    ArrayList<Double> BPoxmf;
    //праві частини
    double[] fFEM;
    double[] fMF;

    public double point;

     ArrayList<JMeshlessInterpolationFunction1D> MeshlessShapeList;

    public MF1SolverLagrangeCoupled(HashMatrix2 Afem, HashMatrix2 Amf, ArrayList<Double> Bfem, ArrayList<Double> Bmf, ArrayList<Double> BPoxfem, ArrayList<Double> BPoxmf, double[] fFEM, double[] fMF) {
        this.Afem = Afem;
        this.Amf = Amf;
        this.Bfem = Bfem;
        this.Bmf = Bmf;
        this.BPoxfem = BPoxfem;
        this.BPoxmf = BPoxmf;
        this.fFEM = fFEM;
        this.fMF = fMF;
    }

    public MF1SolverLagrangeCoupled(HashMatrix2 Afem, HashMatrix2 Amf, ArrayList<Double> Bfem, ArrayList<Double> Bmf, double[] fFEM, double[] fMF, ArrayList<JMeshlessInterpolationFunction1D> MeshlessShapeList) {
        this.Afem = Afem;
        this.Amf = Amf;
        this.Bfem = Bfem;
        this.Bmf = Bmf;
        this.fFEM = fFEM;
        this.fMF = fMF;
        this.MeshlessShapeList = MeshlessShapeList;
    }

    public MF1SolverLagrangeCoupled() {
    }


//    public MF1SolverLagrangeCoupled(HashMatrix Afem, HashMatrix Amf, ArrayList<Double> Bfem, ArrayList<Double> Bmf, double[] fFEM, double[] fMF) {
//        this.Afem = Afem;
//        this.Amf = Amf;
//        this.Bfem = Bfem;
//        this.Bmf = Bmf;
//        this.fFEM = fFEM;
//        this.fMF = fMF;
//    }




    public void generate(){
        int N = Afem.getN()+Amf.getN()+1;
        int n = Amf.getN();

        A = new HashMatrix(HashMatrix.ARRAY_TYPE, N+1);
        b = new double[N+1];
//        безсітковий метод;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                A.setValue(i, j, Amf.getValue(i, j));
            }
            b[i] = fMF[i];
        }
//        МСЕ|
        for (int i = 0; i < Afem.getN(); i++) {
            for (int j = 0; j < Afem.getN(); j++) {
                A.setValue(i+n, j+n, Afem.getValue(i, j));
            }
//            b[i+n] = fFEM[i];
        }
        System.arraycopy(fFEM, 0, b, n, Afem.getN());////////////////////////////////////////////////////////////////////////////////
//        b[Afem.getN()+n] = 0.0;
//        множники Лагранжа
        Bfem.add(0.0);//////////////!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//        for (int i = 0; i < Bmf.size(); i++) {
//            A.setValue(i, N-1, Bmf.get(i));
//            A.setValue(N-1, i, Bmf.get(i));
//        }
//
//        for (int i = 0; i < Bfem.size(); i++) {
//            A.setValue(i+n, N-1, -Bfem.get(i));
//            A.setValue(N-1, i+n, -Bfem.get(i));
//        }
        //point = 0.9;
        for (int i = 0; i < MeshlessShapeList.size(); i++) {
            JMeshlessInterpolationFunction1D jMeshlessInterpolationFunction1D = MeshlessShapeList.get(i);
            double val = jMeshlessInterpolationFunction1D.calculate(point);
            A.setValue(N-1,  Bmf.size() + i, val);
            A.setValue( Bmf.size() + i, N-1, val);
        }
        A.setValue(N-1, Bmf.size()-1, -1);
//        A.setValue(N-1, Bmf.size(), 1);
////        A.setValue(N-1, N-2, -10000);
        A.setValue(Bmf.size()-1, N-1, -1);
//        A.setValue(Bmf.size(), N-1, 1);
////        A.setValue(N-2, N-1, -10000);


//        Правий Діріхле
//        point = 1.0;
//        for (int i = 0; i < MeshlessShapeList.size(); i++) {
//            JMeshlessInterpolationFunction1D jMeshlessInterpolationFunction1D = MeshlessShapeList.get(i);
//            double val = jMeshlessInterpolationFunction1D.calculate(point);
//            if(val<0){continue;}
//            A.setValue(N,  Bmf.size() + i, val);
//            A.setValue( Bmf.size() + i, N, val);
//        }

            point = 1.0;
//            JMeshlessInterpolationFunction1D jMeshlessInterpolationFunction1D = MeshlessShapeList.get(MeshlessShapeList.size()-1);
//            double val = jMeshlessInterpolationFunction1D.calculate(point);
            A.setValue(N,  N-2 , 1);
            A.setValue( N-2 , N, 1);
//
//            jMeshlessInterpolationFunction1D = MeshlessShapeList.get(MeshlessShapeList.size()-2);
//            val = jMeshlessInterpolationFunction1D.calculate(point);
            A.setValue(N,  N-3 , 1);
            A.setValue( N-3 , N, 1);
//            b[N] = 0.01;
    }


    public void generate2(){
        int N = Afem.getN()+Amf.getN()+2;
        int n = Amf.getN();

        A = new HashMatrix(HashMatrix.ARRAY_TYPE, N);
        b = new double[N];
//        безсітковий метод;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                A.setValue(i, j, Amf.getValue(i, j));
            }
            b[i] = fMF[i];
        }
//        МСЕ;
        for (int i = 0; i < Afem.getN(); i++) {
            for (int j = 0; j < Afem.getN(); j++) {
                A.setValue(i+n, j+n, Afem.getValue(i, j));
            }
//            b[i+n] = fFEM[i];
        }
        System.arraycopy(fFEM, 0, b, n, Afem.getN() - 1);///////////////////////////////////////////////////////////////
        b[Afem.getN()+n] = 0.0;
        Bfem.add(0.0);

        for (int i = 0; i < MeshlessShapeList.size(); i++) {
            JMeshlessInterpolationFunction1D jMeshlessInterpolationFunction1D = MeshlessShapeList.get(i);
            double val = jMeshlessInterpolationFunction1D.calculate(point);
            A.setValue(N-1,  Bmf.size() + i, -val);
            A.setValue( Bmf.size() + i, N-1, -val);
        }
//        A.setValue(N-1, Bmf.size()-1, -1);
//        A.setValue(N-1, Bmf.size(), 1);
////        A.setValue(N-1, N-2, -10000);
//        A.setValue(Bmf.size()-1, N-1, -1);


        A.setValue(N-1, Bmf.size()-1, 1);
//        A.setValue(N-1, Bmf.size(), -1);
        A.setValue(Bmf.size(), N-1-1, 1);
//        A.setValue(Bmf.size(), N-1, -1);


        point = 1.0;
        for (int i = 0; i < MeshlessShapeList.size(); i++) {
            JMeshlessInterpolationFunction1D jMeshlessInterpolationFunction1D = MeshlessShapeList.get(i);
            double val = jMeshlessInterpolationFunction1D.calculate(point);
            A.setValue(N-2,  Bmf.size() + i, val);
            A.setValue( Bmf.size() + i, N-2, val);
        }
        b[N-1] = 0.0;
        b[N-2] = dirihletB;
//        A.setValue(N-2, Bmf.size()-1, 1);
//        A.setValue(N-2, Bmf.size(), -1);
//        A.setValue(Bmf.size(), N-2-1, 1);
//        A.setValue(Bmf.size(), N-2, -1);

    }

    public double[] solve(){
//        безсітковий перший фем другий;
        Gauss g = new Gauss(A);
        return  g.calculate(b);
    }

    public CoupledResult solveCoupled(int N/*FEM*/, int nn /*meshFree*/, double divPoint, int Degree, IKernelFunction kernal,double rho){
        ArrayList<Particle> ls1D = new ArrayList<Particle>();
        MeshlessShapeList = new ArrayList<JMeshlessInterpolationFunction1D>();



        double a =divPoint,B=1;

        double h = (B-a)/(double)nn;
        for (int i = 0; i <=nn; i++) {
            ls1D.add(new Particle(a+h*i, 0, i));
            ls1D.get(ls1D.size()-1).setRho(rho*h);
//            MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, a+i*h, 3.1*h, new JGaussian()));
            MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(Degree, a+i*h, rho*h, new JStKernel(rho*h)));
        }

        ArrayList<Double> list1 = new ArrayList<Double>();

        for (int i = 0; i <= nn; i++) {
            list1.add(a+h*i);
        }
        MF1DSolver solver = new MF1DSolver();
            solver.setfString(f);
            solver.setwString(w);
            solver.setqString(q);
            solver.setpString(p);
        solver.setListOfShapeFunction(MeshlessShapeList);
        solver.setNodeList(ls1D);
        solver.setDirihetB(true);
//        dirihletB = 1.0;
        solver.setBVal(0.0);
        solver.setDirihetA(false);
        solver.setAVal(0.0);
        solver.setLeft(divPoint);
        solver.setRight(1.0);


///for FEM !!!!!
        a = 0;
        B = divPoint;
        h = (B-a)/N;
        solver.FEM_H = h;
        solver.generate();
        solver.generatePenalty();
        HashMatrix2 aAa = solver.A;

//        FEM
        ArrayList<Double> list = new ArrayList<Double>();

        for (int i = 0; i <= N; i++) {
            list.add(a+h*i);
        }


        MainMatrixGenerator1 gen = new MainMatrixGenerator1(list);
            gen.setLeft(0.0);
            gen.setRight(divPoint);
            gen.setfString(f);
            gen.setwString(w);
            gen.setqString(q);
            gen.setpString(p);
            gen.setDirihetB(false);
            gen.setBVal(0.0);
            gen.setDirihetA(true);
            gen.setAVal(dirihletA);

            gen.Generate();
            gen.applyBounderiesForCoupling();
        
            MF1SolverLagrangeCoupled cSolverLagrangeCoupled = new MF1SolverLagrangeCoupled(/*solver.getA()*/aAa, gen.getA(),solver.generateBmf(true), gen.generateBFEM(false), solver.getB(),gen.getB(), MeshlessShapeList);
        cSolverLagrangeCoupled.point=divPoint;
        cSolverLagrangeCoupled.generate2();


//        double[] res = cSolverLagrangeCoupled.solve();
        return new CoupledResult(cSolverLagrangeCoupled.solve(), MeshlessShapeList, divPoint, list);

    }

    private ArrayList<Double> generateMidPoints(CoupledResult res){
        ArrayList<Double> list = new ArrayList<Double>();
        for (int i = 0; i < res.getFuncList().size()-1; i++) {
            list.add( (res.getFuncList().get(i).getXi() + res.getFuncList().get(i+1).getXi()) / 2.0);
        }
        
        for (int i = 0; i < res.getMeshlessShapeList().size(); i++) {
//            list.add( (res.getMeshlessShapeList().get(i).getXii() + res.getMeshlessShapeList().get(i+1).getXii()) / 2.0);
        	list.add( (res.getMeshlessShapeList().get(i).getXii() ));
        }

        return list;
    }
    
//    private ArrayList<Double> generateMidPoints(CoupledResult res,int N){
//        ArrayList<Double> list = new ArrayList<Double>();
//        double a = 0.0;
//        double b = 1.0;
//        double h = (b-a)/(double)N;
//        for (int i = 0; i <N; i++) {
//            list.add( a + 3.0*h*i / 2.0);
//        }
//        
//        return list;
//    }

    public double maxNorm(CoupledResult res){
        //1. обчислюєм похідну по середині відрізків-елементів.
                ArrayList<Double> MidPoints = this.generateMidPoints(res);
                ArrayList<Double> derivaties = new ArrayList<Double>();

                for (int i = 0; i < MidPoints.size(); i++) {
                    derivaties.add(res.calculateDerivaty(MidPoints.get(i)));
                }
            //обчислюємо норму
                ArrayList<Double> D = new ArrayList<Double>();
                D.add(Math.abs(res.calculate(MidPoints.get(0)) - res.calculate(MidPoints.get(1))) / Math.abs( MidPoints.get(0) - MidPoints.get(1) ));
                for (int i = 1; i < MidPoints.size()-1; i++) {
                    double d1 = Math.abs(res.calculate(MidPoints.get(i)) - res.calculate(MidPoints.get(i-1))) / Math.abs( MidPoints.get(i) - MidPoints.get(i-1) );
                    double d2 = Math.abs(res.calculate(MidPoints.get(i)) - res.calculate(MidPoints.get(i+1))) / Math.abs( MidPoints.get(i) - MidPoints.get(i+1) );
                    D.add(Math.max(d1, d2));
                    //System.out.println(Math.max(d1, d2));
                }

                D.add(Math.abs(res.calculate(MidPoints.get( MidPoints.size()-1 )) - res.calculate(MidPoints.get(MidPoints.size()-2))) / Math.abs( MidPoints.get(MidPoints.size()-1 ) - MidPoints.get(MidPoints.size()-2) ));
                //System.out.println(Math.abs(res.calculate(MidPoints.get( MidPoints.size()-1 )) - res.calculate(MidPoints.get(MidPoints.size()-2))) / Math.abs( MidPoints.get(MidPoints.size()-1 ) - MidPoints.get(MidPoints.size()-2) ));
                //System.out.println("===========================================");

                double max = -1;
            //обчислюємо Dhk

//                System.out.println("===========================================");
//                for (int k = 0; k < res.getFuncList().size()-1; k++) {
//                    Double d1 = D.get(k);
//                    double hk = res.getFuncList().get(k).gethRight() + res.getFuncList().get(k+1).gethLeft();
//                    double hk = res.getFuncList().get(k+1).getXi() - res.getFuncList().get(k).getXi();
//                    max = Math.max(max, d1*hk);
//                    max = Math.max(max, d1);
//                    D.set(k, d1*hk);
//                    System.out.println(max);
//                }
//                    System.out.println(max);
                for (int k = 0; k < res.getMeshlessShapeList().size()-1; k++) {
                    Double d1 = D.get(k);
                    double hk = res.getMeshlessShapeList().get(k+1).getXii() - res.getMeshlessShapeList().get(k).getXii();
                    max = Math.max(max, d1*hk);
//                    max = Math.max(max, d1);
//                    D.set(k, d1*hk);
                    //System.out.println(max);
                }
//                    System.out.println(max);
//                System.out.println("===========================================");
        return max;
    }


    public double maxError(final CoupledResult res){
    	double max = 0.0;
//    	double[] errr = error2(res);
//    	for (int i = 0; i < errr.length; i++) {
////    		max +=Math.abs(errr[i]);
//			if(Math.abs(errr[i])>max){
//				max = Math.abs(errr[i]);
//			}
//		}
//    	return max/(double)errr.length;
    	max = new GaussQuadrature(50).justGauss(0, 1, new IFunction1D() {
			FEMResult er = erro2(res);
			@Override
			public double calculateDerivaty(double x) {
				// TODO Auto-generated method stub
				return 0;
			}
			
			@Override
			public double calculate(double x) {
				double val = er.calculate(x);
				return val*val;
			}
		});
    	
    	return max;
    }
    
    public double[] error(final CoupledResult femres) {
    	final String pString = p;
    	final String wString = w;
    	final String qString = q;
    	final String fString = f;
    	IFunction1D p = new IFunction1D() {

            public double calculate(double x) {
                try {
                    return PostfiksX.Calculate(pString, x);
                } catch (Exception ex) {
                    Logger.getLogger(MF1DSolver.class.getName()).log(Level.SEVERE, null, ex);
                }
                return 0;
            }

            public double calculateDerivaty(double x) {
                throw new UnsupportedOperationException("Not supported yet.");
            }
        };
        IFunction1D w = new IFunction1D() {

            public double calculate(double x) {
                try {
                    return PostfiksX.Calculate(wString, x);
                } catch (Exception ex) {
                    Logger.getLogger(MF1DSolver.class.getName()).log(Level.SEVERE, null, ex);
                }
                return 0;
            }

            public double calculateDerivaty(double x) {
                throw new UnsupportedOperationException("Not supported yet.");
            }
        };
        IFunction1D q = new IFunction1D() {

            public double calculate(double x) {
                try {
                    return PostfiksX.Calculate(qString, x);
                } catch (Exception ex) {
                    Logger.getLogger(MF1DSolver.class.getName()).log(Level.SEVERE, null, ex);
                }
                return 0;
            }

            public double calculateDerivaty(double x) {
                throw new UnsupportedOperationException("Not supported yet.");
            }
        };
        IFunction1D f = new IFunction1D() {

            public double calculate(double x) {
                try {
                    return PostfiksX.Calculate(fString, x);
                } catch (Exception ex) {
                    Logger.getLogger(MF1DSolver.class.getName()).log(Level.SEVERE, null, ex);
                }
                return 0;
            }

            public double calculateDerivaty(double x) {
                throw new UnsupportedOperationException("Not supported yet.");
            }
        };
    	ArrayList<FEMShapeFunction1D1Pow> FuncList = new ArrayList<FEMShapeFunction1D1Pow>();
    	// for FEM
    	double h1 = 0.01;
    	double h2 = femres.FuncList.get(1).getXi() - femres.FuncList.get(0).getXi();
    	h2/=2.0;
    	FuncList.add(new FEMShapeFunction1D1Pow(0.0, h1, h2));
    	FuncList.add(new FEMShapeFunction1D1Pow(0.0 + h2, h2, h2));
    		
    		
        for (int i = 1; i < femres.FuncList.size()-1; i++) {
        	h1 = femres.FuncList.get(i).getXi() - femres.FuncList.get(i-1).getXi();
        	h1/=2.0;
        	h2 = femres.FuncList.get(i+1).getXi() - femres.FuncList.get(i).getXi();
        	h2/=2.0;
        	FuncList.add(new FEMShapeFunction1D1Pow(femres.FuncList.get(i).getXi(), h1, h2));
        	FuncList.add(new FEMShapeFunction1D1Pow(femres.FuncList.get(i).getXi() + h2, h2, h2));
        }
        
        h1 = femres.FuncList.get(femres.FuncList.size()-1).getXi() - femres.FuncList.get(femres.FuncList.size()-2).getXi();
        h2 = femres.MeshlessShapeList.get(1).getXii() - femres.MeshlessShapeList.get(0).getXii();
    	// for MF
    	
    	h2/=2.0;
    	FuncList.add(new FEMShapeFunction1D1Pow(femres.divPoint, h1, h2));
    	FuncList.add(new FEMShapeFunction1D1Pow(femres.divPoint + h2, h2, h2));
    		
    		
        for (int i = 1; i < femres.FuncList.size()-1; i++) {
        	h1 = femres.MeshlessShapeList.get(i).getXii() - femres.MeshlessShapeList.get(i-1).getXii();
        	h1/=2.0;
        	h2 = femres.MeshlessShapeList.get(i+1).getXii() - femres.MeshlessShapeList.get(i).getXii();
        	h2/=2.0;
        	FuncList.add(new FEMShapeFunction1D1Pow(femres.MeshlessShapeList.get(i).getXii(), h1, h2));
        	FuncList.add(new FEMShapeFunction1D1Pow(femres.MeshlessShapeList.get(i).getXii() + h2, h2, h2));
        }
        
        h1 = femres.MeshlessShapeList.get(femres.MeshlessShapeList.size()-1).getXii() - femres.MeshlessShapeList.get(femres.MeshlessShapeList.size()-2).getXii();
    	h1/=2.0;
    	h2=0.001;
    	FuncList.add(new FEMShapeFunction1D1Pow(femres.MeshlessShapeList.get(femres.MeshlessShapeList.size()-1).getXii(), h1, h2));
    	
        int n = FuncList.size();
		double res[] = null;
		HashMatrix2 A = new HashMatrix2(HashMatrix2.ARRAY_TYPE, n);
		double[] b = new double[n];
		
        
        StatMatrixGenerator gen = new StatMatrixGenerator();
        gen.setLeft(0.0);
        gen.setRight(1.0);

            gen.setF1(FuncList.get(0));
            gen.setF2(FuncList.get(1));
            double[][] tmp = null;
            double[] tmpb = null;

        for (int i = 0; i < b.length-1; i++) {
            gen.setF1(FuncList.get(i));
            gen.setF2(FuncList.get(i+1));
            tmp = gen.getAll(p, q, w);
            tmpb = gen.getFe(f);

            A.setValue(i, i, A.getValue(i, i) + tmp[0][0]);
            A.setValue(i, i+1, A.getValue(i, i+1) + tmp[0][1]);
            A.setValue(i+1, i, A.getValue(i+1, i) + tmp[1][0]);
            A.setValue(i+1, i+1, A.getValue(i+1, i+1) + tmp[1][1]);
            b[i]+=tmpb[0];
            b[i+1]+=tmpb[1];
        }

        b[0] = 0.0;
        A.setValue(0, 1, 0.0);
        A.setValue(0, 0, 1.0);
		b[n-1] = 0.0;
		A.setValue(n-1, n-2, 0.0);
		A.setValue(n-1, n-1, 1.0);
        
        
        double[] tmpU = new double[b.length];
		for (int j = 0; j < tmpU.length; j++) {
			ArrayList<Integer> pos = A.getJRowIndex(j);
			double tmpxz = 0.0;
			for (int i = 0; i < pos.size(); i++) {
				tmpxz += A.getValue(j, pos.get(i))*femres.calculate(FuncList.get(pos.get(i)).getXi());
			}
			b[j]-=tmpxz;
		}

        LUFactorization lu = new LUFactorization(A);
        res = lu.calculate(b);
        XYSeriesCollection col = new XYSeriesCollection();

    		 	XYSeries series = new XYSeries("", true, true);
    	        for (int i = 0; i < FuncList.size(); i++) {
					series.add((double)FuncList.get(i).getXi(),res[i]);
				}
    	        col.addSeries(series);

        NewOkCancelDialog dialog = new NewOkCancelDialog(null, false, col);
        dialog.setLocationRelativeTo(null);
        dialog.setVisible(true);
        
        
		return res;
	}
    
    public double[] error2(final CoupledResult femres) {
    	ArrayList<Double> NodeList = new ArrayList<Double>();
    	int n = 100;
    	double h = (1.0)/n;
        for (int i = 0; i <= n; i++) {
        	NodeList.add(h*i);
        }    	
		double res[] = null;
//		HashMatrix2 A = new HashMatrix2(HashMatrix2.ARRAY_TYPE, n);
//		double[] b = new double[n];
		
        
        MainMatrixGenerator1 gen = new MainMatrixGenerator1(NodeList);
        gen.setLeft(0);
        gen.setRight(1.0);
        gen.setAVal(0.0);
        gen.setfString(f);
        gen.setwString(w);
        gen.setqString(q);
        gen.setpString(p);
        gen.setDirihetB(true);
        gen.setBVal(0.0);
        gen.setDirihetA(true);
        gen.setAVal(0.0);
        
        gen.Generate();
//        gen.applyBounderiesForCoupling();
        
        
        double[] tmpU = new double[gen.b.length];
		for (int j = 0; j < tmpU.length; j++) {
			ArrayList<Integer> pos = gen.A.getJRowIndex(j);
			double tmpxz = 0.0;
			for (int i = 0; i < pos.size(); i++) {
				tmpxz += gen.A.getValue(j, pos.get(i))*femres.calculate(NodeList.get(pos.get(i)));
			}
			gen.b[j]-=tmpxz;
		}

//        LUFactorization lu = new LUFactorization(gen.A);
//        res = lu.calculate(gen.b);
		res = gen.Solve("");
        XYSeriesCollection col = new XYSeriesCollection();

    		 	XYSeries series = new XYSeries("", true, true);
    	        for (int i = 0; i < NodeList.size(); i++) {
					series.add((double)NodeList.get(i),res[i]);
				}
    	        col.addSeries(series);

        NewOkCancelDialog dialog = new NewOkCancelDialog(null, false, col);
        dialog.setLocationRelativeTo(null);
        dialog.setVisible(true);
        
        
		return res;
	}
    
    public FEMResult erro2(final CoupledResult femres) {
    	ArrayList<Double> NodeList = new ArrayList<Double>();
    	int n = 200;
    	double h = (1.0)/n;
        for (int i = 0; i <= n; i++) {
        	NodeList.add(h*i);
        }    	
		double res[] = null;
//		HashMatrix2 A = new HashMatrix2(HashMatrix2.ARRAY_TYPE, n);
//		double[] b = new double[n];
		
        
        nonequidistant.FEMSolver1D.MainMatrixGenerator1 gen = new nonequidistant.FEMSolver1D.MainMatrixGenerator1(NodeList);
        gen.setLeft(0);
        gen.setRight(1.0);
        gen.setAVal(0.0);
        gen.setfString(f);
        gen.setwString(w);
        gen.setqString(q);
        gen.setpString(p);
        gen.setDirihetB(true);
        gen.setBVal(0.0);
        gen.setDirihetA(true);
        gen.setAVal(0.0);
        
        gen.Generate();
//        gen.applyBounderiesForCoupling();
        
        
        double[] tmpU = new double[gen.b.length];
		for (int j = 0; j < tmpU.length; j++) {
			ArrayList<Integer> pos = gen.A.getJRowIndex(j);
			double tmpxz = 0.0;
			for (int i = 0; i < pos.size(); i++) {
				tmpxz += gen.A.getValue(j, pos.get(i))*femres.calculate(NodeList.get(pos.get(i)));
			}
			gen.b[j]-=tmpxz;
		}

//        LUFactorization lu = new LUFactorization(gen.A);
//        res = lu.calculate(gen.b);
		res = gen.Solve("");
        XYSeriesCollection col = new XYSeriesCollection();

    		 	XYSeries series = new XYSeries("", true, true);
    	        for (int i = 0; i < NodeList.size(); i++) {
					series.add((double)NodeList.get(i),res[i]);
				}
    	        col.addSeries(series);

        NewOkCancelDialog dialog = new NewOkCancelDialog(null, false, col);
        dialog.setLocationRelativeTo(null);
        dialog.setVisible(true);
        
        
		return new FEMResult(res, gen.FuncList);
	}
    
    public double Duxotom(double a,double b,double delta,double eps,   int N, int nn, int Degree, IKernelFunction kernal,double rho){
        double x1 , x2;
        double max1 , max2;
        do {
            x1 = (a+b-delta)/2;
            x2 = (a+b+delta)/2;

            CoupledResult res1 = this.solveCoupled(N, nn, x1,  Degree,  kernal, rho);
            CoupledResult res2 = this.solveCoupled(N, nn, x2,  Degree,  kernal, rho);
            max1 = this.maxError(res1);
            max2 = this.maxError(res2);
//            max1 = this.maxNorm(res1);
//            max2 = this.maxNorm(res2);
            if(max1<max2){
                b = x2;
            }else{
                a = x1;
            }

        } while (b-a>eps);
        if(max1<max2){
                return x1;
            }else{
                return x2;
            }
    }
    
    public void cal(double a,double b, double eps,   int N, int nn, int Degree, IKernelFunction kernal,double rho){
    	double x1 = a;
    	while (x1<b) {
    		CoupledResult res1 = this.solveCoupled(N, nn, x1,  Degree,  kernal, rho);
    		
    		
    		
    		XYSeries series1 = new XYSeries("Розвязок", true, true);
    		XYSeries series2 = new XYSeries("Розвязок аналітичний", true, true);
    		double x = 0.0;
//    		System.out.println("f([");
    		while(x<1.001){
    			double tmp = res1.calculate(x);
//    			System.out.println(x+"\t\t"+tmp);
    			series1.add(x,tmp);
    			series2.add(x,(1.0 / (Math.exp(100) - 1) * (1- Math.exp(100*x)) + x));
    			x+=0.01;
    		}
//    		System.out.println("])");
    		XYSeriesCollection col = new XYSeriesCollection(series1);
    		col.addSeries(series2);
    		NewOkCancelDialog dialog = new NewOkCancelDialog(null, false, col);
    		dialog.setLocationRelativeTo(null);
//    		//dialog.setUndecorated(true);
    		dialog.setVisible(true);
    		
    		
    		
            double max1 = this.maxError(res1);
////            System.out.println("===========================================");
            System.out.println(x1+"\t\t" + max1);
            x1+=0.01;
		}
    }

    public double GoldenDivision(double a,double b, double eps,   int N, int nn, int Degree, IKernelFunction kernal,double rho){
        double x1 , x2;
        double max1 , max2;
        double x = a + (3-Math.sqrt(5.0))/2*(b-a);
        int numm = 0;
        CoupledResult resprev1 = this.solveCoupled(N, nn, 0.5,  Degree,  kernal, rho);
        CoupledResult resprev2 =resprev1;
        do {
        	numm++;
            x1 = a + (3-Math.sqrt(5.0))/2*(b-a);
            x2 = a + (Math.sqrt(5.0)-1)/2*(b-a);

            final CoupledResult res1 = this.solveCoupled(N, nn, x1,  Degree,  kernal, rho);
            final CoupledResult res2 = this.solveCoupled(N, nn, x2,  Degree,  kernal, rho);
//            max1 = this.maxNorm(res1);
//            max2 = this.maxNorm(res2);
            
//            max1 = this.maxError(res1);
//            max2 = this.maxError(res2);
            final CoupledResult res11 = resprev1;
            GaussQuadrature gau = new GaussQuadrature(100);
        	max1 = gau.Gauss(0.0, 1.0, 0.1, new IFunction1D() {
        		
                public double calculate(double x) {
                	double result1 = res11.calculate(x);
                	double result2 = res1.calculate(x);
                	return (result1-result2)*(result1-result2);
                }

                public double calculateDerivaty(double x) {
                    throw new UnsupportedOperationException("Not supported yet.");
                }
            });
        	final CoupledResult res22 = resprev2;
        	max2 = gau.Gauss(0.0, 1.0, 0.1, new IFunction1D() {
        		
                public double calculate(double x) {
                	double result1 = res22.calculate(x);
                	double result2 = res2.calculate(x);
                	return (result1-result2)*(result1-result2);
                }

                public double calculateDerivaty(double x) {
                    throw new UnsupportedOperationException("Not supported yet.");
                }
            });
        	resprev1 = res1;
        	resprev2 = res2;
            System.out.println("===========================================");
            System.out.println(max1+"\t\t" + x1);
            System.out.println(max2+"\t\t" + x2);

//            double xxx = 0.0;
//            System.out.println("===========================================");
//            while (xxx<1.0000001) {
//                System.out.println(xxx+"\t\t"+res1.calculate(xxx));
//                xxx+=0.01;
//            }
//            System.out.println(numm+"===========================================");

            if(max1<max2){
                b = x2;
                if(Math.abs(x-x1)>Math.abs(x-x2)){
                    x = x2;
                }else{
                    x = x1;
                }

            }else{
                a = x1;
                if(Math.abs(x-x1)>Math.abs(x-x2)){
                    x = x2;
                }else{
                    x = x1;
                }
            }

        } while (b-a>eps);
        System.out.println("к-сть ітерацій = "+numm);
        return x;
    }
    
    public double GoldenDivision2(double a,double b, double eps,   int N, int nn, int Degree, IKernelFunction kernal,double rho){
        double x1 , x2;
        double max1 , max2;
        double x = a + (3-Math.sqrt(5.0))/2*(b-a);
        int numm = 0;
        do {
        	numm++;
            x1 = a + (3-Math.sqrt(5.0))/2*(b-a);
            x2 = a + (Math.sqrt(5.0)-1)/2*(b-a);

            final CoupledResult res1 = this.solveCoupled(N, nn, x1,  Degree,  kernal, rho);
            final CoupledResult res2 = this.solveCoupled(N, nn, x2,  Degree,  kernal, rho);
            max1 = this.maxError(res1);
            max2 = this.maxError(res2);
            
//            max1 = this.maxError(res1);
//            max2 = this.maxError(res2);
            System.out.println("===========================================");
            System.out.println(max1+"\t\t" + x1);
            System.out.println(max2+"\t\t" + x2);

//            double xxx = 0.0;
//            System.out.println("===========================================");
//            while (xxx<1.0000001) {
//                System.out.println(xxx+"\t\t"+res1.calculate(xxx));
//                xxx+=0.01;
//            }
//            System.out.println(numm+"===========================================");

            if(max1<max2){
                b = x2;
                if(Math.abs(x-x1)>Math.abs(x-x2)){
                    x = x2;
                }else{
                    x = x1;
                }

            }else{
                a = x1;
                if(Math.abs(x-x1)>Math.abs(x-x2)){
                    x = x2;
                }else{
                    x = x1;
                }
            }

        } while (b-a>eps);
        System.out.println("к-сть ітерацій = "+numm);
        return x;
    }
    
    public optimizationResult getOptimal(double eps,double divPoint, int nFEM,int nEFG, int Degree, IKernelFunction kernal,double rho){
    	boolean key = true;
    	int Num = 0;
    	int nFEMtmp = 0;
    	int nEFGtmp = 0;
    	double max1 = 10.0;
    	double maxPrev = 0.0;
    	GaussQuadrature gau = new GaussQuadrature(100);
    	while(true){
    		Num++;
    		
    		
    		if(key){
    			key=false;
    			nFEMtmp = nFEM;
    			nEFGtmp = nEFG;    			
    		}
    		
   			// Збиваємо тимчасові дані по н для перерахунку опитмальних ...
//    			nFEMtmp = nFEM;
//    			nEFGtmp = nEFG; 
    		//1 оптимізуємо по МСЕ
    			CoupledResult res1 = this.solveCoupled(nFEMtmp, nEFGtmp, divPoint,  Degree,  kernal, rho);
//                max1 = this.maxNorm(res1);
                double max2 = 0.0;
                
                do {
                	nFEMtmp+=5;
                	final CoupledResult res2 = this.solveCoupled(nFEMtmp, nEFGtmp, divPoint,  Degree,  kernal, rho);
                	final CoupledResult resprev = res1;
                	max2 = gau.Gauss(0.0, 1.0, 0.1, new IFunction1D() {
                		
                        public double calculate(double x) {
                        	double result1 = resprev.calculate(x);
                        	double result2 = res2.calculate(x);
                        	return (result1-result2)*(result1-result2);
                        }

                        public double calculateDerivaty(double x) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }
                    });
                	res1=res2;
//                	max2 = this.maxNorm(res1);
                	System.out.println("max1\t"+max1+"\tmax2\t"+max2+"\tdelta\t"+Math.abs(max2-max1));
                	System.out.println("fem\t"+nFEMtmp);
                	if((max1<max2)){nFEMtmp-=5; break;}
                	if(Math.abs(max1-max2)<eps){
                		max1 = max2; 
                		nFEMtmp-=5;
                		break;
                	}else{
                		max1 = max2;
                	}
				} while (true);
    		//2 оптимізуємо по EFG
                do {
                	nEFGtmp+=5;
                	final CoupledResult res2 = this.solveCoupled(nFEMtmp, nEFGtmp, divPoint,  Degree,  kernal, rho);
                	final CoupledResult resprev = res1;
                	max2 = gau.Gauss(0.0, 1.0, 0.1, new IFunction1D() {
                		
                        public double calculate(double x) {
                        	double result1 = resprev.calculate(x);
                        	double result2 = res2.calculate(x);
                        	return (result1-result2)*(result1-result2);
                        }

                        public double calculateDerivaty(double x) {
                            throw new UnsupportedOperationException("Not supported yet.");
                        }
                    });
                	res1=res2;
//                	max2 = this.maxNorm(res1);
                	System.out.println("EFG\t"+nEFGtmp);
                	System.out.println("max1\t"+max1+"max2\t"+max2);
                	if((max1<max2)){nEFGtmp-=5; break;}
                	if(Math.abs(max1-max2)<eps){
                		max1 = max2;
                		nEFGtmp-=5;
                		break;
                	}else{
                		max1 = max2;
                	}

    		
				} while (true);
                System.out.println("==============================================");
                System.out.println(Num+"\tFEM = "+nFEMtmp+"\tEFG = "+nEFGtmp+"\tdivPoint = "+divPoint+"\tmax = "+max1);
            //0 оптимізуємо по точці розбиття
//    			divPoint = GoldenDivision(divPoint-0.1, divPoint+0.1, eps, nFEMtmp, nEFGtmp, Degree, kernal, rho);
//    		divPoint = GoldenDivision(0.8, 0.99, eps, nFEMtmp, nEFGtmp, Degree, kernal, rho);
    		divPoint = GoldenDivision(0.5, 0.99, 0.01, nFEMtmp, nEFGtmp, Degree, kernal, rho);
    		if(Math.abs(max1-maxPrev)<eps*Math.abs(maxPrev)){
            		break;
            	}else{
            		maxPrev = max1;                		
            	}
    	}
    	return new optimizationResult(nFEMtmp, nEFGtmp, divPoint);
    }
    
    public optimizationResult getOptimal2(double eps,double divPoint, int nFEM,int nEFG, int Degree, IKernelFunction kernal,double rho){
    	boolean key = true;
    	int Num = 0;
    	int nFEMtmp = 0;
    	int nEFGtmp = 0;
    	double max1 = 10.0;
    	double maxPrev = 0.0;
    	while(true){
    		Num++;
    		
    		
    		if(key){
    			key=false;
    			nFEMtmp = nFEM;
    			nEFGtmp = nEFG;    			
    		}
    		
   			// Збиваємо тимчасові дані по н для перерахунку опитмальних ...
//    			nFEMtmp = nFEM;
//    			nEFGtmp = nEFG; 
    		//1 оптимізуємо по МСЕ
    			CoupledResult res1 = this.solveCoupled(nFEMtmp, nEFGtmp, divPoint,  Degree,  kernal, rho);
                max1 = this.maxError(res1);
                double max2 = 0.0;
                
                do {
                	nFEMtmp+=5;
                	res1 = this.solveCoupled(nFEMtmp, nEFGtmp, divPoint,  Degree,  kernal, rho);
                	max2 = this.maxError(res1);
                	System.out.println("max1\t"+max1+"\tmax2\t"+max2+"\tdelta\t"+Math.abs(max2-max1));
                	System.out.println("fem\t"+nFEMtmp);
                	if((max1<max2)){nFEMtmp-=5; break;}
                	if(Math.abs(max1-max2)<eps){
//                		max1 = max2; 
                		nFEMtmp-=5;
                		break;
                	}else{
                		max1 = max2;
                	}
				} while (true);
    		//2 оптимізуємо по EFG
                do {
                	nEFGtmp+=5;
                	res1 = this.solveCoupled(nFEMtmp, nEFGtmp, divPoint,  Degree,  kernal, rho);
                	max2 = this.maxError(res1);
                	System.out.println("EFG\t"+nEFGtmp);
                	System.out.println("max1\t"+max1+"\tmax2"+max2);
                	if((max1<max2)){nEFGtmp-=5; break;}
                	if(Math.abs(max1-max2)<eps){
//                		max1 = max2;
                		nEFGtmp-=5;
                		break;
                	}else{
                		max1 = max2;
                	}

    		
				} while (true);
                System.out.println("==============================================");
                System.out.println(Num+"\tFEM = "+nFEMtmp+"\tEFG = "+nEFGtmp+"\tdivPoint = "+divPoint+"\tmax = "+max1);
              //0 оптимізуємо по точці розбиття
//        		divPoint = GoldenDivision2(0.8, 0.99, 0.01, nFEMtmp, nEFGtmp, Degree, kernal, rho);
        		divPoint = Duxotom(0.8, 0.999, eps, 0.01, nFEMtmp, nEFGtmp, Degree, kernal, rho);
//    			divPoint = GoldenDivision(divPoint-0.1, divPoint+0.1, eps, nFEMtmp, nEFGtmp, Degree, kernal, rho);
//    		divPoint = GoldenDivision(0.8, 0.99, eps, nFEMtmp, nEFGtmp, Degree, kernal, rho);
    		
                if(Math.abs(divPoint-maxPrev)<eps){
            		break;
            	}else{
            		maxPrev = divPoint;                		
            	}
    	}
    	return new optimizationResult(nFEMtmp, nEFGtmp, divPoint);
    }

    public CoupledResult autoCoupling(double Start, double End, int N, int nn,int nnn, int Degree, IKernelFunction kernal,double rho){
        ArrayList<Double> divPoints = new ArrayList<Double>();
        CoupledResult res = null;
//        int nnn = 10;
        double h = (End-Start)/(double)nnn;
        for (int i = 0; i <= nnn; i++) {
            divPoints.add(Start+h*i);
        }
        double min = 10000000;
        for (int i = 0; i < divPoints.size(); i++) {
            Double divPointLocal = divPoints.get(i);
//            this.solveCoupled(N, nn, divPoint);
//            this.maxNorm(this.solveCoupled(N, nn, divPoint));
            final CoupledResult tmpres = this.solveCoupled(N, nn, divPointLocal,  Degree,  new JStKernel(rho), rho);
            double tmp = this.maxNorm(tmpres);
            if(tmp<min){
                min = tmp;
                res = tmpres;
                divPoint = divPointLocal;
            }
//            System.out.println("point : "+divPointLocal+"\t norma : "+tmp+"");
//            System.out.println(+divPointLocal+"\t"+tmp+"");
            GaussQuadrature gau = new GaussQuadrature(100);
            System.out.println("\t\t"+gau.Gauss(0.0, 1.0, 0.1, new IFunction1D() {

                public double calculate(double x) {
                	double result = tmpres.calculate(x);
                	return result*result;
                }

                public double calculateDerivaty(double x) {
                    throw new UnsupportedOperationException("Not supported yet.");
                }
            }));

        }
        return res;
    }
    
    public CoupledResult auto(double Start, double End, int N, int nn,int nnn, int Degree, IKernelFunction kernal,double rho){
        ArrayList<Double> divPoints = new ArrayList<Double>();
        CoupledResult res = null;
//        int nnn = 10;
        double h = (End-Start)/(double)nnn;
        for (int i = 0; i <= nnn; i++) {
            divPoints.add(Start+h*i);
        }
        double min = 10000000;
        for (int i = 0; i < divPoints.size(); i++) {
            Double divPointLocal = divPoints.get(i);
//            this.solveCoupled(N, nn, divPoint);
//            this.maxNorm(this.solveCoupled(N, nn, divPoint));
            final CoupledResult tmpres = this.solveCoupled(N, nn, divPointLocal,  Degree,  new JStKernel(rho), rho);
            final CoupledResult tmpresBetter = this.solveCoupled(2*N, 2*nn, divPointLocal,  Degree,  new JStKernel(rho), rho);
//            double tmp = this.maxNorm(tmpres);
            GaussQuadrature gau = new GaussQuadrature(100);
            double tmp = gau.Gauss(0.0, 1.0, 100, new IFunction1D() {
            	
            	public double calculate(double x) {
            		double result = tmpres.calculate(x);
            		double result2 = tmpresBetter.calculate(x);
            		return (result-result2)*(result-result2);
            	}

            	public double calculateDerivaty(double x) {
            		throw new UnsupportedOperationException("Not supported yet.");
            	}
            });

//            double tmp1 = gau.Gauss(0.0, 1.0, 100, new IFunction1D() {
//
//                public double calculate(double x) {
//                	double result = tmpres.calculate(x);
//                	return result*result;
//                }
//
//                public double calculateDerivaty(double x) {
//                    throw new UnsupportedOperationException("Not supported yet.");
//                }
//            });
//
//            double tmp2 = gau.Gauss(0.0, 1.0, 100, new IFunction1D() {
//
//                public double calculate(double x) {
//                	double result = tmpresBetter.calculate(x);
//                	return result*result;
//                }
//
//                public double calculateDerivaty(double x) {
//                    throw new UnsupportedOperationException("Not supported yet.");
//                }
//            });
//            System.out.println(tmp1+"\t"+tmp2+"\t"+Math.abs(tmp1-tmp2));
//            if(Math.abs(tmp1-tmp2)<min){
            System.out.println(tmp);
            if(Math.abs(tmp)<min){
                min = Math.abs(tmp);
                res = tmpres;
                divPoint = divPointLocal;
            }
//            System.out.println("point : "+divPointLocal+"\t norma : "+tmp+"");
//            System.out.println(+divPointLocal+"\t"+tmp+"");
//            GaussQuadrature gau = new GaussQuadrature(100);
//            System.out.println("\t\t"+gau.Gauss(0.0, 1.0, 0.1, new IFunction1D() {
//
//                public double calculate(double x) {
//                	double result = tmpres.calculate(x);
//                	return result*result;
//                }
//
//                public double calculateDerivaty(double x) {
//                    throw new UnsupportedOperationException("Not supported yet.");
//                }
//            }));

        }
        return res;
    }

    static double divPoint = 0.929937500;//0.955470937;//0.97;//0.84875;//0.868;0.8383089635127199;//

    public static void main(String[] args){
        int nn = 5; // meshFree
        int N = 5; // FEM
        f = "100*sin(3.14*x)";
        q = "0";
        w = "100";
        p = "1";
        dirihletB = 1;
        dirihletA = 0;
        
//        double divPoint = 0.7475;//0.9464670950385429;
        MF1SolverLagrangeCoupled s = new MF1SolverLagrangeCoupled();
//        s.autoCoupling(0.5, 0.995, N, nn,10,1,new JGaussian(1),1.1);
//        s.auto(0.8, 0.995, N, nn,20,1,new JGaussian(1),2.1);
        
        System.out.println(divPoint);

        ArrayList<Particle> ls1D = new ArrayList<Particle>();
        final ArrayList<JMeshlessInterpolationFunction1D> MeshlessShapeList = new ArrayList<JMeshlessInterpolationFunction1D>();

//        System.out.println("========================================");
//        System.out.println("=========дихотомія======================");
//        divPoint = s.Duxotom(0.5, 0.999, 0.001, 0.01, N, nn, 1, new JGaussian(), 1.1);
//        System.out.println(divPoint);
//
//        System.out.println("========================================");
//        System.out.println("=========золотий поділ==================");
//        divPoint = s.GoldenDivision(0.5, 0.999, 1E-5, N, nn, 1, new JGaussian(1.1), 1.1);
//        System.out.println("Результат = "+divPoint);
//        System.out.println("========================================");
      System.out.println("=========cal==================");
//        s.cal(0.88, 0.999, 1E-5, N, nn, 1, new JGaussian(1.1), 2.1);
        
        System.out.println("========================================");
        System.out.println("=========optimizationResult==================");
//        optimizationResult or = s.getOptimal(0.001, divPoint, N, nn, 1, new JGaussian(1), 2.1);  0.00000001
//        optimizationResult or = s.getOptimal2(0.0000000001, divPoint, N, nn, 1, new JGaussian(1), 2.1);
//        N = or.getN(); 						  
//        nn = or.getM();
//        divPoint = or.getDivPoint();
//        System.out.println("FEM = "+N+"\tEFG = "+nn+"\tdivPoint = "+divPoint);

        double a =divPoint,b=1;
        double rho = 2.1;
        double h = (b-a)/(double)nn;
        for (int i = 0; i <=nn; i++) {
            ls1D.add(new Particle(a+h*i, 0, i));
            ls1D.get(ls1D.size()-1).setRho(rho*h);
            //MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, a+i*h, 1.33*h, new JExponential(0.3)));
            MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, a+i*h, rho*h, new JStKernel(rho*h)));
            //MeshlessShapeList.add(new JMeshlessInterpolationFunction1D(1, a+i*h, 1.33*h, new JStKernel()));
            //System.out.print(ls1D.get(i).getX()+"\t");
        }
        


        ArrayList<Double> list1 = new ArrayList<Double>();

        for (int i = 0; i <= nn; i++) {
            list1.add(a+h*i); 
            
        }
        
        for (int i = 0; i < MeshlessShapeList.size(); i++) {
        	MeshlessShapeList.get(i).setNodeList(ls1D);
		}
        
        XYSeriesCollection col = new XYSeriesCollection();
        for (Iterator<JMeshlessInterpolationFunction1D> it = MeshlessShapeList.iterator(); it.hasNext();) {
    		JMeshlessInterpolationFunction1D jm = (JMeshlessInterpolationFunction1D) it.next();
    		double x1 = a;
    		 	XYSeries series = new XYSeries("", true, true);
    	        while(x1<=1.0001){
    	            double tmp1 = jm.calculate(x1);
    	            series.add(x1,tmp1);
    	            x1+=0.002;
    	        }
    	        col.addSeries(series);
    	}
        NewOkCancelDialog dialog = new NewOkCancelDialog(null, false, col);
        dialog.setLocationRelativeTo(null);
        dialog.setVisible(true);

//        MainMatrixGenerator1 solver = new MainMatrixGenerator1(list1);
//            solver.setLeft(0.5);
//            solver.setRight(1.0);
//            solver.setfString("100");
//            solver.setwString("100");
//            solver.setqString("0");
//            solver.setpString("1");
//            solver.setDirihetB(true);
//            solver.setBVal(0.0);
//            solver.setDirihetA(false);
//            solver.setAVal(0.0);
//
//
//            solver.Generate();
//            solver.applyBounderiesForCoupling();

//        MF1DSolver solver = new MF1DSolver();
//            solver.setfString("100");
//            solver.setwString("100");
//            solver.setqString("0");
//            solver.setpString("1");
//        solver.setListOfShapeFunction(MeshlessShapeList);
//        solver.setNodeList(ls1D);
//        solver.setDirihetB(true);
//        solver.setBVal(0.0);
//        solver.setDirihetA(false);
//        solver.setAVal(0.0);
//        solver.setLeft(divPoint);
//        solver.setRight(1.0);
        
        MF1DSolver solver = new MF1DSolver();
        solver.setfString(f);
        solver.setwString(w);
        solver.setqString(q);
        solver.setpString(p);
        solver.setAVal(0.0);
        solver.setBVal(1.0);
        solver.setDirihetA(false);
        solver.setDirihetB(true);
//        dirihletB = 1.0;
        solver.setRight(b);
        solver.setLeft(a);

        solver.setListOfShapeFunction(MeshlessShapeList);
        solver.setNodeList(ls1D);


        long startTime = System.currentTimeMillis();


        solver.generate();


        ///for FEM !!!!!
        a = 0;
        b = divPoint;
        
        h = (b-a)/N;



        solver.FEM_H = h;
//
//
        solver.generate();
//
//
//        solver.generatePenalty();
        HashMatrix2 aAa = solver.A;
//         aAa = solver.generateLagrangeMulti();
        



//        FEM
        ArrayList<Double> list = new ArrayList<Double>();
        
        for (int i = 0; i <= N; i++) {
            list.add(a+h*i);
        }


        MainMatrixGenerator1 gen = new MainMatrixGenerator1(list);
            gen.setLeft(0.0);
            gen.setRight(divPoint);
            gen.setfString(f);
            gen.setwString(w);
            gen.setqString(q);
            gen.setpString(p);
            gen.setDirihetB(false);
            gen.setBVal(0.0);
            gen.setDirihetA(true);
            gen.setAVal(dirihletA);

            gen.Generate();
            gen.applyBounderiesForCoupling();
        //gen.Solve("");


        //MF1SolverLagrangeCoupled cSolverLagrangeCoupled = new MF1SolverLagrangeCoupled(gen.getA(), solver.getA(),gen.generateBFEM(false), solver.generateBmf(true), gen.getB(), solver.getB());
            //MF1SolverLagrangeCoupled cSolverLagrangeCoupled = new MF1SolverLagrangeCoupled(gen.getA(), solver.getA(),gen.generateBFEM(false), solver.generateBFEM(true), gen.getB(), solver.getB());
        //MF1SolverLagrangeCoupled cSolverLagrangeCoupled = new MF1SolverLagrangeCoupled(solver.getA(), gen.getA(),solver.generateBmf(true), gen.generateBFEM(false), solver.getB(), gen.getB());
//            MF1SolverLagrangeCoupled cSolverLagrangeCoupled = new MF1SolverLagrangeCoupled(solver.getA(), gen.getA(),solver.generateBmf(true), gen.generateBFEM(false),solver.generateBmf(true), gen.generateBFEMPox(false), solver.getB(), gen.getB());
//        MF1SolverLagrangeCoupled cSolverLagrangeCoupled = new MF1SolverLagrangeCoupled(gen.getA(), solver.getA(),gen.generateBFEM(false), solver.generateBmf(true), gen.getB(), solver.getB(),MeshlessShapeList);

            MF1SolverLagrangeCoupled cSolverLagrangeCoupled = new MF1SolverLagrangeCoupled(/*solver.getA()*/aAa, gen.getA(),solver.generateBmf(true), gen.generateBFEM(false), solver.getB(),gen.getB(), MeshlessShapeList);
        cSolverLagrangeCoupled.point=divPoint;
//        cSolverLagrangeCoupled.generate();
            cSolverLagrangeCoupled.generate2();

//        for (int i = 0; i < cSolverLagrangeCoupled.b.length; i++) {
//            for (int j = 0; j < cSolverLagrangeCoupled.b.length; j++) {
//                System.out.print(cSolverLagrangeCoupled.A.getValue(i, j)+"\t\t\t");
//            }
//            System.out.println("\t\t\t"+cSolverLagrangeCoupled.b[i]);
//        }
        
        System.out.println();
        System.out.println("=========================================================");



//        double[] res = solver.solveLagrangeMulti();
        final double[] res = cSolverLagrangeCoupled.solve();

//        System.out.print("===================================================");
//        for (int i = 0; i < res.length; i++) {
//            System.out.print(res[i]+"\n");
//        }
//        System.out.print("===================================================");

        XYSeries series1 = new XYSeries("Розвязок", true, true);
        XYSeries series2 = new XYSeries("Розвязок2", true, true);
        CoupledResult xs = new CoupledResult(res, MeshlessShapeList, divPoint, list);
        h = 0.001;
        for ( double x = 0; x <= 1+h; x+=h) {
        	
        	double nu = 0.01;
        	double A = 1;
        	double aux = Math.PI*(A*A+nu*nu*Math.PI*Math.PI);
            double e = Math.exp(A/nu);
            double c1 = (-aux+A*(e+1))/(aux*(e-1));
            double c2 = (aux-2*A)/(aux*(e-1));
            double r = c1 + c2*Math.exp(A*x/nu) + nu*Math.PI*(Math.sin(Math.PI*x)-A*Math.cos(Math.PI*x)/(nu*Math.PI))/aux;
        	
        	double tmp = -(1.0 / (Math.exp(100) - 1) * (1- Math.exp(100*x)) );
        	tmp = r;
//        	tmp = (1.0 / (100*Math.exp(100)) * (1- Math.exp(100*x)) + x);
        	series2.add(x,tmp);
//            series2.add(x,xs.calculate(x));
        }

        
//        int NuM = list.size();
        for (int i = 0; i < list.size(); i++) {
            System.out.print(list.get(i)+"\t"+res[i]+"\n");
            series1.add((double)list.get(i),(double)res[i]);
        }
//        for (int i = 0; i < list1.size(); i++) {
//            System.out.print(list1.get(i)+"\t"+res[i+NuM]+"\n");
//            series1.add((double)list1.get(i),(double)res[i+NuM]);
//        }

        b=1.0;
        h = (b-divPoint)/(1000.0);
        
        for ( double x = divPoint+h; x <= b+h; x+=h) {
            double tmp = 0.0;
            for (int j = 0; j <MeshlessShapeList.size(); j++) {
                tmp+= MeshlessShapeList.get(j).calculate(x)*res[j+N+1];
            }
            System.out.print(x+"\t"+tmp+"\n");
            series1.add(x,tmp);
        }

//        while(x<=1.001){
//            double tmp = 0.0;
//            for (int j = 0; j <MeshlessShapeList.size(); j++) {
//                tmp+= MeshlessShapeList.get(j).calculate(x)*res[j+N];
//            }
//            System.out.print(x+"\t"+tmp+"\n");
//            series1.add(x,tmp);
//            x+=h;
//        }
        System.out.println("Результат = "+divPoint);
        final CoupledResult cr = new CoupledResult(res, MeshlessShapeList, divPoint, list);
        GaussQuadrature gau = new GaussQuadrature(100);
        System.out.println("FEM = "+N+"\tEFG = "+nn+"\tdivPoint = "+divPoint);
        System.out.println("L2\t\t"+gau.Gauss(0.0, 1.0, 100, new IFunction1D() {

            public double calculate(double x) {
            	double result = cr.calculate(x);
            	return result*result;
            }

            public double calculateDerivaty(double x) {
                throw new UnsupportedOperationException("Not supported yet.");
            }
        }));
        
        System.out.println("L2 err\t\t"+gau.Gauss(0.0, 1.0, 100, new IFunction1D() {

            public double calculate(double x) {
            	double result = cr.calculate(x);
            	double nu = 0.01;
            	double A = 1;
            	double aux = Math.PI*(A*A+nu*nu*Math.PI*Math.PI);
                double e = Math.exp(A/nu);
                double c1 = (-aux+A*(e+1))/(aux*(e-1));
                double c2 = (aux-2*A)/(aux*(e-1));
                double r = c1 + c2*Math.exp(A*x/nu) + nu*Math.PI*(Math.sin(Math.PI*x)-A*Math.cos(Math.PI*x)/(nu*Math.PI))/aux;
            	result-=r;
//            	result -= -(1.0 / (Math.exp(100) - 1) * (1- Math.exp(100*x)) );
            	
            	return result*result;
            }

            public double calculateDerivaty(double x) {
                throw new UnsupportedOperationException("Not supported yet.");
            }
        }));
        
        
            
//        0.8383089635127199
//        L2		0.3136698991974013
        System.out.println();
//        System.out.println("L2 real\t\t"+gau.Gauss(0.0, 1.0, 0.1, new IFunction1D() {
//
//            public double calculate(double x) {
//            	double result = (1.0 / (Math.exp(100) - 1) * (1- Math.exp(100*x)) + x );
//            	return result*result;
//            }
//
//            public double calculateDerivaty(double x) {
//                throw new UnsupportedOperationException("Not supported yet.");
//            }
//        }));
        System.out.println("=========================================================");
        col = new XYSeriesCollection(series1);
        col.addSeries(series2);
        dialog = new NewOkCancelDialog(null, false, col);
        dialog.setLocationRelativeTo(null);
        //dialog.setUndecorated(true);
        dialog.setVisible(true);

//        double x = 0.0;
//        XYSeries series = new XYSeries("розвязок", true, true);
//        x = 0.0;
//        while(x<=1.0001){
//            double tmp = 0.0;
//            tmp = MeshlessShapeList.get(0).calculate(x);
////            for (int j = 0; j <MeshlessShapeList.size(); j++) {
////                tmp+= MeshlessShapeList.get(j).calculate(x)*res[j];
////            }
//            System.out.print(x+"\t"+tmp+"\n");
//            series.add(x,tmp);
//            x+=0.01;
//        }
//        NewOkCancelDialog dialog1 = new NewOkCancelDialog(null, false, new XYSeriesCollection(series));
//            dialog1.setLocationRelativeTo(null);
//            //dialog.setUndecorated(true);
//            dialog1.setVisible(true);

    }

    public static JMeshlessInterpolationFunctionLiu jm;
    public static JMeshlessInterpolationFunction1D jm1D;
    
}
